/*   2013/01/29   */
// Run III Data Compression tool
// read raw data file and read/write compressed data file
#include <iostream>
#include <fstream>
#include <sstream>
#include <string>
#include <time.h>
#include <TSystem.h>
using namespace std;

#include "TF1.h"
#include "TH1.h"
#include "TCanvas.h"
#include "TROOT.h"
#include "TMath.h"
#include "TStyle.h"
#include "TGraph.h"
#include "TApplication.h"

// A : atmic mass
// Mt : target mass
// Mwimp : wimp mass
// sigmawn : Xsec for wimp and nucleon

double c = 3.0e8; // [m/s]
double ck = 3.0e5; // [km/s]
double rho0 = 0.3; //[GeV/c2/cm3] WIMP energy density
double Mwimp = 100; // [GeV/c2] WIMP mass
double Mn = 0.9315; // [GeV/c2] average nucleon mass
double v0 = 220;    // [km/s] Sun's velocity from Galactic center
double vesc = 544;  // [km/s] Galactic escape velocity
//double vesc = 600;  // [km/s] Galactic escape velocity
double vear = 232;  // [km/s] Average earth velocity from Galactic center
double NA = 6.0221415e23;  // Avogadro's number
double pi = TMath::Pi();

double c1 = 0.751;      //factor change for each months
double c2 = 0.561;      //factor change for each months

double Mr(double M1, double M2); // reduced mass
double Rr(double M1, double M2); // reduced mass
double F(double ER, double A, double mT);


double dRdE(double ER ,double A , double sigmawn){
  // argon 

  double Mt = A * Mn; // GeV/c2
  double rm = Rr(Mwimp,Mt);
  double vmin = sqrt( 2.0*ER / rm / Mwimp )*c/1000.;

  //  printf("%f %f %f %f %f %f %f\n",ER*1.0e6,Mwimp,vesc,vmin,vear,vesc-vear,vesc+vear);

  if (vesc+vear < vmin) return 0.;  
  
  double sigmaA = pow(A,2) * pow( Mr(Mwimp,Mt)/Mr(Mwimp,Mn),2 ) * sigmawn;
  double R0 = (361/Mwimp/Mt) * (sigmaA) * (rho0/0.3) * (v0/220); 
  double kfk0 =  TMath::Erf(vesc/v0) - 2/sqrt(pi)*vesc/v0*exp(-pow(vesc/v0,2) ) ;
  double E0 = Mwimp*pow((v0/ck),2)/2.;         // GeV*(km/s)^2   E(100GeV)=mc^2
  
  if (vmin > 0 && vmin < vesc-vear) {  

    double Term1 = sqrt(pi)*v0/(4*vear) * ( TMath::Erf((vmin+vear)/v0) - TMath::Erf((vmin-vear)/v0) );
    double Term2 = exp( -pow(vesc/v0,2) );
    double drate = R0*kfk0/E0/rm * (Term1 - Term2) * 1.0e-6;
    return drate ;
  }

  if (vmin > vesc - vear  && vmin < vesc + vear) {  

    double Term1 = sqrt(pi)*v0/(4*vear) * ( TMath::Erf(vesc/v0) - TMath::Erf((vmin-vear)/v0) );
    double Term2 = (vesc+vear-vmin)/2.0/vear*exp(-pow(vesc/v0,2) );
    double drate = R0*kfk0/E0/rm * (Term1 - Term2) * 1.0e-6;
    return drate ;

  }
}

double dRdEold(double ER ,double A , double sigmawn){
  // argon 
  double Mt = A * 0.932; // GeV/c2
  double sigmaA = pow(A,2) * pow( Mr(Mwimp,Mt)/Mr(Mwimp,Mn),2 ) * sigmawn;
  double vmin = sqrt( Mt*ER / (2*Mr(Mwimp,Mt)) );
  double Q = sqrt(2*Mt*ER)/c; // keV s/m

  double k0 = pow(pi*pow(v0,2),3/2);
  double kf = k0 * ( TMath::Erf(vesc/v0) - 2/sqrt(pi)*vesc/v0*exp( -pow(vesc/v0,2) ) );
  //  double R0 = 2 *NA * rho0/sqrt(pi)/A/Mwimp * sigmaA * v0;
  double R0 = (361/Mwimp/Mt) * (sigmaA) * (rho0/0.3) * (v0/220);   // GeV , pb , GeV/cm3 , km/s  ## Nishimura ph.D ##

  double E0 = Mwimp*pow((v0/ck),2)/2.;         // GeV*(km/s)^2   E(100GeV)=mc^2
  double rm = 4*Mwimp*Mt/pow(Mwimp+Mt,2); 

  //  double Term1 = sqrt(pi)*v0/(4*vear) * ( TMath::Erf((vmin+vear)/v0) - TMath::Erf((vmin-vear)/v0) );
  double Term1 = c1 * exp(-c2 * ER / E0 / rm);  // year average
  double Term2 = exp( -pow(vesc/v0,2) );
 
  //***1
  //  double drate =  R0/E0/rm * exp(-ER/E0/rm) *1e-6;  // vesc = inf , vear = 0
  //***2
  double drate =  R0/E0/rm * Term1 *1e-6;  // vesc = inf , vear = vear
  //***3
  //double drate = k0/kf * R0/E0/rm * ( Term1 - Term2 ) *1e-6; //  vesc = vesc , vear = vear 

    //  return  R0/E0/rm * exp(-ER/E0/rm) *1e-6;
  return drate;//100kg*100day
}

// reduce mass
double Mr(double M1, double M2){
  return M1*M2/(M1+M2);
}

// reduce mass
double Rr(double M1, double M2){
  return 4.0*M1*M2/(M1+M2)/(M1+M2);
}


//form factor  Helm form factor (CDMS thesis)           

double F(double ER, double A, double mT){ // ER -> recoil E , mT -> target mass (keV)
  int option=0;
  
  double Q = sqrt(2*mT*ER)/c; // keV s/m                                              

  double aT = 0.52                    *pow(10,-15); // m                              
  double cT = (1.23*pow(A,1/3.) - 0.6)*pow(10,-15); // m                              
  double sT = 0.9                     *pow(10,-15); // m                              

  double rT = sqrt( pow(cT,2)+7/3.*pow(pi*aT,2)-5*pow(sT,2) ); // [m] effective nuclear radius  

  double hbar = 6.582e-19; // keVs                                                    
  double xq = Q*rT/hbar; // dimentionless                                             
  
  double formfac = 3*(sin(xq) - xq*cos(xq))/pow(xq,3) * exp( -pow(Q*sT/hbar,2)/2. );       
  if(option==1) formfac = 3*sqrt(pi/(2*xq))*TMath::BesselJ1(xq)/xq * exp( -pow(Q*sT/hbar,2)/2. );  // use TMath BesselJ1 if you want 

 return formfac;
}

//-----------------                                                                   
// test form factor                                                                   
//-----------------                                                                   
double FormFac(double *x, double *par){
  double ER1 = x[0];
  double A1 = par[0];
  double mT1 = par[1];
  return pow(F(ER1,A1,mT1),2);
}

// for plot  -> dRdE * FormFactor

double dForm(double *x, double *par){
  double er = x[0];
  double A  = par[0];
  double sigmawn  = par[1];
  double mT = par[2];
  
  double er2=er * 1e-6;
  double drF = dRdE(er2,A,sigmawn) * pow(F(er,A,mT),2);
  //  double drF = dRdE(er2,A,sigmawn);
  //  return 100.*100.*drF;
  return drF;
}

double dFormold(double *x, double *par){
  double er = x[0];
  double A  = par[0];
  double sigmawn  = par[1];
  double mT = par[2];
  
  double er2=er * 1e-6;
  double drF = dRdEold(er2,A,sigmawn) * pow(F(er,A,mT),2);
  //  double drF = dRdE(er2,A,sigmawn);
  //  return 100.*100.*drF;
  return drF;
}

void drdeF(double Mw){

  Mwimp = Mw;

  TF1 *f2 = new TF1("f2", dForm, 0, 120, 3);
  f2->SetParameters(40,1e-6,40*0.932e6);    // A->Ar(40), sigma , mT =A*0.932(keV)
  TF1 *f2old = new TF1("f2", dFormold, 0, 120, 3);
  f2old->SetParameters(40,1e-6,40*0.932e6);    // A->Ar(40), sigma , mT =A*0.932(keV)

  TCanvas *can = new TCanvas("can", "can");
  //  can->SetLogx();
  can->SetLogy();
  //f1->SetMinimum(1e-7)
  //f1->SetTitle("dRdE");
  //f1->SetLineColor(2);
  //f1->Draw("");
  f2->GetXaxis()->SetTitle("Recoil Energy [keV]");
  f2->GetYaxis()->SetTitle("Count [kev/kg/day]");


  f2->SetMinimum(1e-7);
  f2->SetMaximum(1e-1);
  f2->SetTitle("dRdE");
  f2->SetLineColor(4);
  f2old->SetLineColor(2);
  f2->Draw("");
  f2old->Draw("same");

}


void lim(){

  TF1 *f2 = new TF1("f2", dForm, 0, 120, 3);
  f2->SetParameters(40,1.0,40*0.932e6);    // A->Ar(40), sigma , mT =A*0.932(keV)
  TH1F *h1 = new TH1F("h1","",1000,0.,1000.);
  double gx[10][1000];
  double gy[10][1000];
  float  mass[10]={100.,10000.,100.,100.,100.,10000.};
  float  ecut[10]={60.,60.,20.,5.,5.,5.};
  int ng[6] = {0};
  for (int i = 0;i<1000;i++){
    Mwimp = (float) (i+1.0);
    for (int j = 0;j<6;j++){
      float ncut = f2->Integral(ecut[j],200.);
      if (ncut > 0){
	double l90 = 1.0/(1.0*ncut)*2.3;
	printf("%i %f %f %e\n",j,Mwimp,ncut,l90);
	gx[j][ng[j]]=Mwimp;
	gy[j][ng[j]]=l90 * 1.0e-36 / mass[j];
	ng[j]++;
      }
    }
  }

  h1->GetXaxis()->SetRange(5,1000);
  h1->SetMinimum(1.0e-45);
  h1->SetMaximum(1.0e-39);
  h1->SetXTitle("WIMP Mass (GeV/c)");
  h1->SetYTitle("WIMP-Nucleon Cross Section (cm^{2})");
  h1->Draw();
  TGraph *gg[6];
  for (int j=0;j<4;j++){
    gg[j] = new TGraph(ng[j],gx[j],gy[j]);
    if (j==0) {gg[j]->SetLineColor(1);}
    if (j==1) {gg[j]->SetLineColor(4);}
    if (j==2) gg[j]->SetLineColor(1);
    if (j==3) gg[j]->SetLineColor(1);
    if (j==4) gg[j]->SetLineColor(1);
    if (j==5) gg[j]->SetLineColor(4);
    gg[j]->Draw("Lsame");
  }
}

